home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Power Programmierung
/
Power-Programmierung (Tewi)(1994).iso
/
magazine
/
progjour
/
1988
/
06
/
ccode
/
fltbench.c
< prev
Wrap
Text File
|
1988-09-01
|
39KB
|
1,465 lines
#ifdef __HIGHC__
#pragma wtrsvd1(16); /* New Intel setting for Weitek registers. */
#endif
/* ------------------------------------------------------------------
Floating-point benchmarks.
-------------------------------------------------------------------*/
#include <stdio.h>
#include <math.h>
#include "bench.h"
#define TRUE 1
#define FALSE 0
/*
main - Main driver module for benchmarks
*/
null_init(){}
main(L_argc,L_argv)
int L_argc; char ** L_argv;
{
extern Float();
extern flt();
extern savage();
extern autodesk(), autodesk_init(), INTRIG_autodesk();
extern whetd(),whets();
long Wtime;
argc = L_argc; argv = L_argv;
printf("\n");
printf("Floating-point benchmarks from BYTE Magazine (July 1987, Page 101),\n");
printf("and other standard benchmarks -- version 1.1, corrected Whetstones.\n");
printf("\n");
common_header(); /* Dobench knows the format of this. */
DoBench("F", "Float", 10000L, Float, null_init, TRUE);
DoBench("f", "Flt", 256000L, flt, null_init, TRUE);
DoBench("S", "Savage", 25000L, savage, null_init, TRUE);
DoBench("A", "Autodesk",1000L, autodesk, autodesk_init, TRUE);
DoBench("AI", "Autodesk (INTRIG)",1000L, INTRIG_autodesk, autodesk_init, TRUE);
Wtime = DoBench("WS", "Whetstone (single)", 100L, whets, null_init, FALSE);
if (Wtime) printf(" (%d KWhets/sec).\n",(int)(10000 / (Wtime/100.0)));
Wtime = DoBench("WD", "Whetstone (double)", 100L, whetd, null_init, FALSE);
if (Wtime) printf(" (%d KWhets/sec).\n",(int)(10000 / (Wtime/100.0)));
common_trailer();
}
/*------------------------------------------------------------------*/
/* simple benchmark for testing floating point speed of c libraries
does repeated multiplications and divisions in a loop that is
large enough to make the looping time insignificant */
#define CONST1 3.141597E0
#define CONST2 1.7839032E4
#define COUNT 10000
Float()
{
double a, b, c;
int i;
a = CONST1;
b = CONST2;
for (i = 0; i < COUNT; ++i)
{
c = a * b;
c = c / a;
c = a * b;
c = c / a;
c = a * b;
c = c / a;
c = a * b;
c = c / a;
c = a * b;
c = c / a;
c = a * b;
c = c / a;
c = a * b;
c = c / a;
}
/* printf ("Done\n"); */
}
#undef COUNT
/*------------------------------------------------------------------*/
/***************/
/* Flt.c */
/***************/
flt() {
long i,j,loops; double a,b,c;
loops = 256000;
for( i=1 ; i<=loops ; i++ ) {
j = loops - i;
a = (double)i;
b = (double)j;
c = b / a;
a = b - c;
b = c * a;
c = b + a;
}
a = c + b;
}
/*------------------------------------------------------------------*/
/*
** savage.c -- floating point speed and accuracy test. C version
** derived from BASIC version which appeared in Dr. Dobb's Journal,
** Sep. 1983, pp. 120-122.
*/
#define ILOOP 25000
extern double tan(), atan(), exp(), log(), sqrt();
savage()
{
int i;
double a;
/* printf("start\n"); */
a = 1.0;
for (i = 1; i <= (ILOOP - 1); i++)
a = tan(atan(exp(log(sqrt(a*a))))) + 1.0;
/* printf("a = %20.14e\n", a); */
/* printf("done\n"); */
}
/*------------------------------------------------------------------*/
/**************************************************************************
* *
* Whetstone benchmark in C. This program is a translation of the *
* original Algol version in "A Synthetic Benchmark" by H.J. Curnow *
* and B.A. Wichman in Computer Journal, Vol 19 #1, February 1976. *
* *
* Used to test compiler efficiency, optimization, and double *
* precision floating-point performance. This version is specific *
* to the Turbo-Amiga and Amiga but it can be easily adapted to *
* other systems by replacing the clock() routine with your own. *
* *
**************************************************************************/
#define ITERATIONS 10 /* 1 Million Whetstone instructions */
/* #define POUT */ /* Leave as is for 'Official' result. */
/* #define MTEST */ /* define for Module timing results */
/* only. Leave as is for 'Official' */
/* result. */
static double D_e1[4];
static double D_t, D_t1, D_t2;
static long j, k, l;
whetd() {
long i;
long n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11;
long m, loops;
double D_x, D_y, D_z, D_x1, D_x2, D_x3, D_x4;
/************************/
/* initialize constants */
/************************/
D_t = 0.499975;
D_t1 = 0.500250;
D_t2 = 2.0;
/***********************/
/* Set Module Weights. */
/***********************/
m = 10; /* m = 10 is used to obtain better timing */
loops = m * ITERATIONS; /* accuracy only. Slow systems should use */
n1 = 0 * loops; /* m = 1. */
n2 = 12 * loops;
n3 = 14 * loops;
n4 = 345 * loops;
n5 = 0 * loops;
n6 = 210 * loops;
n7 = 32 * loops;
n8 = 899 * loops;
n9 = 616 * loops;
n10 = 0 * loops;
n11 = 93 * loops;
/*********************************/
/* MODULE 1: simple identifiers */
/*********************************/
D_x1 = 1.0;
D_x2 = -1.0;
D_x3 = -1.0;
D_x4 = -1.0;
if( n1 > 0 )
{
for(i = 1; i <= n1; i++)
{
D_x1 = ( D_x1 + D_x2 + D_x3 - D_x4 ) * D_t;
D_x2 = ( D_x1 + D_x2 - D_x3 - D_x4 ) * D_t;
D_x3 = ( D_x1 - D_x2 + D_x3 + D_x4 ) * D_t;
D_x4 = (-D_x1 + D_x2 + D_x3 + D_x4 ) * D_t;
}
}
/*****************************/
/* MODULE 2: Array Elements */
/*****************************/
D_e1[0] = 1.0; /* Start at element 0 in C, vice 1 in Fortran */
D_e1[1] = -1.0;
D_e1[2] = -1.0;
D_e1[3] = -1.0;
if( n2 > 0 )
{
for (i = 1; i <= n2; i++)
{
D_e1[0] = ( D_e1[0] + D_e1[1] + D_e1[2] - D_e1[3] ) * D_t;
D_e1[1] = ( D_e1[0] + D_e1[1] - D_e1[2] + D_e1[3] ) * D_t;
D_e1[2] = ( D_e1[0] - D_e1[1] + D_e1[2] + D_e1[3] ) * D_t;
D_e1[3] = (-D_e1[0] + D_e1[1] + D_e1[2] + D_e1[3] ) * D_t;
}
}
/*********************************/
/* MODULE 3: Array as Parameter */
/*********************************/
if( n3 > 0 )
{
for (i = 1; i <= n3; i++)
{
D_pa(D_e1);
}
}
/********************************/
/* MODULE 4: Conditional Jumps */
/********************************/
j = 1;
if( n4 > 0 )
{
for (i = 1; i <= n4; i++)
{
if (j == 1)
j = 2;
else
j = 3;
if (j > 2)
j = 0;
else
j = 1;
if (j < 1 )
j = 1;
else
j = 0;
}
}
/**********************/
/* MODULE 5: Omitted */
/**********************/
/*********************************/
/* MODULE 6: Integer Arithmetic */
/*********************************/
j = 1;
k = 2;
l = 3;
if( n6 > 0 )
{
for (i = 1; i <= n6; i++)
{
j = j * (k - j) * (l -k);
k = l * k - (l - j) * k;
l = (l - k) * (k + j);
D_e1[l - 2] = j + k + l; /* Remember we started at D_e1[0]. */
D_e1[k - 2] = j * k * l; /* l-2 in C, vice l-1 in Fortran */
}
}
/**************************************/
/* MODULE 7: Trigonometric Functions */
/**************************************/
D_x = 0.5;
D_y = 0.5;
if( n7 > 0 )
{
for(i = 1; i <= n7; i++)
{
D_x = D_t * atan(D_t2*sin(D_x)*cos(D_x)/(cos(D_x+D_y)+cos(D_x-D_y)-1.0));
D_y = D_t * atan(D_t2*sin(D_y)*cos(D_y)/(cos(D_x+D_y)+cos(D_x-D_y)-1.0));
}
}
/******************************/
/* MODULE 8: Procedure Calls */
/******************************/
D_x = 1.0;
D_y = 1.0;
D_z = 1.0;
if( n8 > 0 )
{
for (i = 1; i <= n8; i++)
{
D_p3(D_x, D_y, &D_z);
}
}
/*******************************/
/* MODULE 9: Array References */
/*******************************/
j = 1;
k = 2;
l = 3;
D_e1[0] = 1.0;
D_e1[1] = 2.0;
D_e1[2] = 3.0;
if( n9 > 0 )
{
for(i = 1; i <= n9; i++)
{
D_p0();
}
}
/**********************************/
/* MODULE 10: Integer Arithmetic */
/**********************************/
j = 2;
k = 3;
if( n10 > 0 )
{
for(i = 1; i <= n10; i++)
{
j = j + k;
k = j + k;
j = k - j;
k = k - j - j;
}
}
/**********************************/
/* MODULE 11: Standard Functions */
/**********************************/
D_x = 0.75;
if( n11 > 0 )
{
for(i = 1; i <= n11; i++)
{
D_x = sqrt( exp( log(D_x) / D_t1) );
}
}
/**************************/
/* End of Whetstone Tests */
/**************************/
#if 0 /* Done elsewhere now. */
stoptime = clock();
benchtime = stoptime - starttime - nulltime;
D_x1 = (double)benchtime/100.0;
printf(" Benchtime(sec) = %lf\n",D_x1);
D_x2 = 100.0 * (double)loops / D_x1;
KWhets = (long)D_x2;
printf(" KWhets/sec = %ld\n\n",KWhets);
#endif
}
/*******************/
/* Subroutine pa() */
/*******************/
D_pa(e) /* Exactly as in the Algol 60 version, but we */
/* could do away with that 'goto'. */
double e[4];
{
int j;
j = 0;
lab:
e[0] = ( e[0] + e[1] + e[2] - e[3] ) * D_t;
e[1] = ( e[0] + e[1] - e[2] + e[3] ) * D_t;
e[2] = ( e[0] - e[1] + e[2] + e[3] ) * D_t;
e[3] = ( -e[0] + e[1] + e[2] + e[3] ) / D_t2;
j ++;
if (j < 6)
goto lab;
}
/************************/
/* Subroutine p3(x,y,z) */
/************************/
D_p3(x, y, z)
double x, y, *z;
{
x = D_t * (x + y);
y = D_t * (x + y);
*z = (x + y) /D_t2;
}
/*******************/
/* Subroutine p0() */
/*******************/
D_p0()
{
D_e1[j] = D_e1[k];
D_e1[k] = D_e1[l];
D_e1[l] = D_e1[j];
}
static float S_e1[4];
static float S_t, S_t1, S_t2;
/* A true single-precision whestone can only be achieved by using
ANSI prototypes. Without them, doubles get passed to p3 and converted
to float upon entry to p3.
If your compiler doesn't support prototypes, you could instead
pass pointers to x and y, for a slight decrease in efficiency.
*/
#ifndef NO_PROTOTYPES
S_p3(float,float,float*);
#endif
whets() {
long i;
long n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11;
long m, loops;
float S_x, S_y, S_z, S_x1, S_x2, S_x3, S_x4;
/************************/
/* initialize constants */
/************************/
S_t = 0.499975;
S_t1 = 0.500250;
S_t2 = 2.0;
/***********************/
/* Set Module Weights. */
/***********************/
m = 10; /* m = 10 is used to obtain better timing */
loops = m * ITERATIONS; /* accuracy only. Slow systems should use */
n1 = 0 * loops; /* m = 1. */
n2 = 12 * loops;
n3 = 14 * loops;
n4 = 345 * loops;
n5 = 0 * loops;
n6 = 210 * loops;
n7 = 32 * loops;
n8 = 899 * loops;
n9 = 616 * loops;
n10 = 0 * loops;
n11 = 93 * loops;
/*********************************/
/* MODULE 1: simple identifiers */
/*********************************/
S_x1 = 1.0;
S_x2 = -1.0;
S_x3 = -1.0;
S_x4 = -1.0;
if( n1 > 0 )
{
for(i = 1; i <= n1; i++)
{
S_x1 = ( S_x1 + S_x2 + S_x3 - S_x4 ) * S_t;
S_x2 = ( S_x1 + S_x2 - S_x3 - S_x4 ) * S_t;
S_x3 = ( S_x1 - S_x2 + S_x3 + S_x4 ) * S_t;
S_x4 = (-S_x1 + S_x2 + S_x3 + S_x4 ) * S_t;
}
}
/*****************************/
/* MODULE 2: Array Elements */
/*****************************/
S_e1[0] = 1.0; /* Start at element 0 in C, vice 1 in Fortran */
S_e1[1] = -1.0;
S_e1[2] = -1.0;
S_e1[3] = -1.0;
if( n2 > 0 )
{
for (i = 1; i <= n2; i++)
{
S_e1[0] = ( S_e1[0] + S_e1[1] + S_e1[2] - S_e1[3] ) * S_t;
S_e1[1] = ( S_e1[0] + S_e1[1] - S_e1[2] + S_e1[3] ) * S_t;
S_e1[2] = ( S_e1[0] - S_e1[1] + S_e1[2] + S_e1[3] ) * S_t;
S_e1[3] = (-S_e1[0] + S_e1[1] + S_e1[2] + S_e1[3] ) * S_t;
}
}
/*********************************/
/* MODULE 3: Array as Parameter */
/*********************************/
if( n3 > 0 )
{
for (i = 1; i <= n3; i++)
{
S_pa(S_e1);
}
}
/********************************/
/* MODULE 4: Conditional Jumps */
/********************************/
j = 1;
if( n4 > 0 )
{
for (i = 1; i <= n4; i++)
{
if (j == 1)
j = 2;
else
j = 3;
if (j > 2)
j = 0;
else
j = 1;
if (j < 1 )
j = 1;
else
j = 0;
}
}
/**********************/
/* MODULE 5: Omitted */
/**********************/
/*********************************/
/* MODULE 6: Integer Arithmetic */
/*********************************/
j = 1;
k = 2;
l = 3;
if( n6 > 0 )
{
for (i = 1; i <= n6; i++)
{
j = j * (k - j) * (l -k);
k = l * k - (l - j) * k;
l = (l - k) * (k + j);
S_e1[l - 2] = j + k + l; /* Remember we started at S_e1[0]. */
S_e1[k - 2] = j * k * l; /* l-2 in C, vice l-1 in Fortran */
}
}
/**************************************/
/* MODULE 7: Trigonometric Functions */
/**************************************/
S_x = 0.5;
S_y = 0.5;
if( n7 > 0 )
{
for(i = 1; i <= n7; i++)
{
S_x = S_t * atan(S_t2*sin(S_x)*cos(S_x)/(cos(S_x+S_y)+cos(S_x-S_y)-1));
S_y = S_t * atan(S_t2*sin(S_y)*cos(S_y)/(cos(S_x+S_y)+cos(S_x-S_y)-1));
}
}
/******************************/
/* MODULE 8: Procedure Calls */
/******************************/
S_x = 1.0;
S_y = 1.0;
S_z = 1.0;
if( n8 > 0 )
{
for (i = 1; i <= n8; i++)
{
S_p3(S_x, S_y, &S_z);
}
}
/*******************************/
/* MODULE 9: Array References */
/*******************************/
j = 1;
k = 2;
l = 3;
S_e1[0] = 1.0;
S_e1[1] = 2.0;
S_e1[2] = 3.0;
if( n9 > 0 )
{
for(i = 1; i <= n9; i++)
{
S_p0();
}
}
/**********************************/
/* MODULE 10: Integer Arithmetic */
/**********************************/
j = 2;
k = 3;
if( n10 > 0 )
{
for(i = 1; i <= n10; i++)
{
j = j + k;
k = j + k;
j = k - j;
k = k - j - j;
}
}
/**********************************/
/* MODULE 11: Standard Functions */
/**********************************/
S_x = 0.75;
if( n11 > 0 )
{
for(i = 1; i <= n11; i++)
{
S_x = sqrt( exp( log(S_x) / S_t1) );
}
}
/**************************/
/* End of Whetstone Tests */
/**************************/
#if 0 /* Done elsewhere now. */
stoptime = clock();
benchtime = stoptime - starttime - nulltime;
S_x1 = (double)benchtime/100.0;
printf(" Benchtime(sec) = %lf\n",S_x1);
S_x2 = 100.0 * (double)loops / S_x1;
KWhets = (long)S_x2;
printf(" KWhets/sec = %ld\n\n",KWhets);
#endif
}
/*******************/
/* Subroutine pa() */
/*******************/
S_pa(e) /* Exactly as in the Algol 60 version, but we */
/* could do away with that 'goto'. */
float e[4];
{
int j;
j = 0;
lab:
e[0] = ( e[0] + e[1] + e[2] - e[3] ) * S_t;
e[1] = ( e[0] + e[1] - e[2] + e[3] ) * S_t;
e[2] = ( e[0] - e[1] + e[2] + e[3] ) * S_t;
e[3] = ( -e[0] + e[1] + e[2] + e[3] ) / S_t2;
j ++;
if (j < 6)
goto lab;
}
/************************/
/* Subroutine p3(x,y,z) */
/************************/
S_p3(x, y, z)
float x, y, *z;
{
x = S_t * (x + y);
y = S_t * (x + y);
*z = (x + y) /S_t2;
}
/*******************/
/* Subroutine p0() */
/*******************/
S_p0()
{
S_e1[j] = S_e1[k];
S_e1[k] = S_e1[l];
S_e1[l] = S_e1[j];
}
/*-- End Of Whetstone C Source Code -----------------*/
/*------------------------------------------------------------------*/
/* Autodesk benchmark. */
/*------------------------------------------------------------------*/
/*
This is a complete optical design raytracing algorithm,
stripped of its user interface and recast into portable C. It
not only determines execution speed on an extremely floating
point (including trig function) intensive real-world
application, it checks accuracy on an algorithm that is
exquisitely sensitive to errors. The performance of this
program is typically far more sensitive to changes in the
efficiency of the trigonometric library routines than the
average floating point program.
The benchmark may be compiled in two modes. If the symbol
INTRIG is defined, built-in trigonometric and square root
routines will be used for all calculations. Timings made with
INTRIG defined reflect the machine's basic floating point
performance for the arithmetic operators. If INTRIG is not
defined, the system library <math.h> functions are used. Results
with INTRIG not defined reflect the system's library performance
and/or floating point hardware support for trig functions and
square root. Results with INTRIG defined are a good guide to
general floating point performance, while results with INTRIG
undefined indicate the performance of an application which is
math function intensive.
Special note regarding errors in accuracy: this program has
generated numbers identical to the last digit it formats and
checks on the following machines, floating point architectures,
and languages:
IBM PC / XT / AT Lattice C IEEE 64 bit, 80 bit temporaries
High C same, in line 80x87 code
BASICA "Double precision"
Quick BASIC IEEE double precision, software routines
Sun 3 C IEEE 64 bit, 80 bit temporaries,
in-line 68881 code, in-line FPA code.
MicroVAX II C Vax "G" format floating point
Macintosh Plus MPW C SANE floating point, IEEE 64 bit format
implemented in ROM.
Inaccuracies reported by this program should be taken VERY
SERIOUSLY INDEED, as the program has been demonstrated to be
invariant under changes in floating point format, as long as
the format is a recognised double precision format. If you
encounter errors, please remember that they are just as likely
to be in the floating point editing library or the
trigonometric libraries as in the low level operator code.
The benchmark assumes that results are basically reliable, and
only tests the last result computed against the reference. If
you're running on a suspect system you can compile this program
with ACCURACY defined. This will generate a version which
executes as an infinite loop, performing the ray trace and
checking the results on every pass. All incorrect results will
be reported.
Representative timings are given below. All have been normalised
as if run for 1000 iterations.
Time in seconds Computer, Compiler, and notes
Normal INTRIG
3466.00 4031.00 Commodore 128, 2 Mhz 8510 with software floating
point. Abacus Software/Data-Becker Super-C 128,
version 3.00, run in fast (2 Mhz) mode. Note:
the results generated by this system differed
from the reference results in the 8th to 10th
decimal place.
3290.00 IBM PC/AT 6 Mhz, Microsoft/IBM BASICA version A3.00.
Run with the "/d" switch, software floating point.
2131.50 IBM PC/AT 6 Mhz, Lattice C version 2.14, small model.
This version of Lattice compiles subroutine
calls which either do software floating point
or use the 80x87. The machine on which I ran
this had an 80287, but the results were so bad
I wonder if it was being used.
1598.00 Macintosh Plus, MPW C, SANE Software floating point.
404.00 IBM PC/AT 6 Mhz, Microsoft QuickBASIC version 2.0.
Software floating point.
165.15 IBM PC/AT 6 Mhz, Metaware High C version 1.3, small
model. This was compiled to call subroutines for
floating point, and the machine contained an 80287
which was used by the subroutines.
143.20 Macintosh II, MPW C, SANE calls. I was unable to
determine whether SANE was using the 68881 chip or
not.
121.80 Sun 3/160 16 Mhz, Sun C. Compiled with -fsoft switch
which executes floating point in software.
78.78 IBM RT PC (Model 6150). IBM AIX 1.0 C compiler
with -O switch.
66.36 206.35 IBM PC/AT 6Mhz, Metaware High C version 1.3, small
model. This was compiled with in-line code for the
80287 math coprocessor. Trig functions still call
library routines.
17.18 Apollo DN-3000, 12 Mhz 68020 with 68881, compiled
with in-line code for the 68881 coprocessor.
According to Apollo, the library routines are chosen
at runtime based on coprocessor presence. Since the
coprocessor was present, the library is supposed to
use in-line floating point code.
15.55 27.56 VAXstation II GPX. Compiled and executed under
VAX/VMS C.
15.14 37.93 Macintosh II, Unix system V. Green Hills 68020
Unix compiler with in-line code for the 68881
coprocessor (-O -ZI switches).
12.69 Sun 3/160 16 Mhz, Sun C. Compiled with -fswitch,
which calls a subroutine to select the fastest
floating point processor. This was using the 68881.
11.74 26.73 Compaq Deskpro 386, 16 Mhz 80386 with 16 Mhz 80387.
Metaware High C version 1.3, compiled with in-line
for the math coprocessor (but not optimised for the
80386/80387). Trig functions still call library
routines.
8.43 30.49 Sun 3/160 16 Mhz, Sun C. Compiled with -f68881,
generating in-line MC68881 instructions. Trig
functions still call library routines.
6.29 25.17 Sun 3/260 25 Mhz, Sun C. Compiled with -f68881,
generating in-line MC68881 instructions. Trig
functions still call library routines.
6.00 23.00 Compaq Deskpro 386, 16 Mhz 80386 with 16 Mhz 80387.
Metaware High C version 1.3, compiled with in-line
for the math coprocessor (optimised for the
80386/80387).
5.9 23.00 Compaq Deskpro 386, 16 Mhz 80386 with 16 Mhz 80387.
Metaware High C version 1.4, compiled with in-line
for the math coprocessor (optimised for the
80386/80387).
5.25 --- Intel 386/20, 16 Mhz 80386 with 16 Mhz 80387.
Metaware High C version 1.4, compiled with in-line
for the math coprocessor (optimised for the
80386/80387).
*/
#define cot(x) (1.0 / tan(x))
#define INTRIG_cot(x) (1.0 / INTRIG_tan(x))
#define max_surfaces 10
/* Local variables */
static char tbfr[132];
static short current_surfaces;
static short paraxial;
static double clear_aperture;
static double aberr_lspher;
static double aberr_osc;
static double aberr_lchrom;
static double max_lspher;
static double max_osc;
static double max_lchrom;
static double radius_of_curvature;
static double object_distance;
static double ray_height;
static double axis_slope_angle;
static double from_index;
static double to_index;
static double spectral_line[9];
static double s[max_surfaces][5];
static double od_sa[2][2];
static char outarr[8][80]; /* Computed output of program goes here */
int itercount; /* The iteration counter for the main loop
in the program is made global so that
the compiler should not be allowed to
optimise out the loop over the ray
tracing code. */
/* The test case used in this program is the design for a 4 inch
achromatic telescope objective used as the example in Wyld's
classic work on ray tracing by hand, given in Amateur Telescope
Making, Volume 3. */
static double testcase[4][4] = {
{27.05, 1.5137, 63.6, 0.52},
{-16.68, 1, 0, 0.138},
{-16.68, 1.6164, 36.7, 0.38},
{-78.1, 1, 0, 0}
};
/* Internal trig functions (used only if INTRIG is defined). These
standard functions may be enabled to obtain timings that reflect
the machine's floating point performance rather than the speed of
its trig function evaluation. */
#define INTRIG_fabs(x) ((x < 0.0) ? -x : x)
#define pic 3.1415926535897932
/* Commonly used constants */
static double pi = pic,
twopi =pic * 2.0,
piover4 = pic / 4.0,
fouroverpi = 4.0 / pic,
piover2 = pic / 2.0;
/* Coefficients for ATAN evaluation */
static double atanc[] = {
0.0,
0.4636476090008061165,
0.7853981633974483094,
0.98279372324732906714,
1.1071487177940905022,
1.1902899496825317322,
1.2490457723982544262,
1.2924966677897852673,
1.3258176636680324644
};
/* aint(x) Return integer part of number. Truncates towards 0 */
double aint(x)
double x;
{
long l;
/* Note that this routine cannot handle the full floating point
number range. This function should be in the machine-dependent
floating point library! */
l = x;
if ((int)(-0.5) != 0 && l < 0 )
l++;
x = l;
return x;
}
/* sin(x) Return sine, x in radians */
double INTRIG_sin(x)
double x;
{
int sign;
double y, r, z;
x = (((sign= (x < 0.0)) != 0) ? -x: x);
if (x > twopi)
x -= (aint(x / twopi) * twopi);
if (x > pi) {
x -= pi;
sign = !sign;
}
if (x > piover2)
x = pi - x;
if (x < piover4) {
y = x * fouroverpi;
z = y * y;
r = y * (((((((-0.202253129293E-13 * z + 0.69481520350522E-11) * z -
0.17572474176170806E-8) * z + 0.313361688917325348E-6) * z -
0.365762041821464001E-4) * z + 0.249039457019271628E-2) * z -
0.0807455121882807815) * z + 0.785398163397448310);
} else {
y = (piover2 - x) * fouroverpi;
z = y * y;
r = ((((((-0.38577620372E-12 * z + 0.11500497024263E-9) * z -
0.2461136382637005E-7) * z + 0.359086044588581953E-5) * z -
0.325991886926687550E-3) * z + 0.0158543442438154109) * z -
0.308425137534042452) * z + 1.0;
}
return sign ? -r : r;
}
/* cos(x) Return cosine, x in radians, by identity */
double INTRIG_cos(x)
double x;
{
x = (x < 0.0) ? -x : x;
if (x > twopi) /* Do range reduction here to limit */
x = x - (aint(x / twopi) * twopi); /* roundoff on add of PI/2 */
return INTRIG_sin(x + piover2);
}
/* tan(x) Return tangent, x in radians, by identity */
static double INTRIG_tan(x)
double x;
{
return INTRIG_sin(x) / INTRIG_cos(x);
}
/* sqrt(x) Return square root. Initial guess, then Newton-
Raphson refinement */
double INTRIG_sqrt(x)
double x;
{
double c, cl, y;
int n;
if (x == 0.0)
return 0.0;
if (x < 0.0) {
fprintf(stderr,
"\nGood work! You tried to take the square root of %g",
x);
fprintf(stderr,
"\nunfortunately, that is too complex for me to handle.\n");
exit(1);
}
y = (0.154116 + 1.893872 * x) / (1.0 + 1.047988 * x);
c = (y - x / y) / 2.0;
cl = 0.0;
for (n = 50; c != cl && n--;) {
y = y - c;
cl = c;
c = (y - x / y) / 2.0;
}
return y;
}
/* atan(x) Return arctangent in radians,
range -pi/2 to pi/2 */
double INTRIG_atan(x)
double x;
{
int sign, l, y;
double a, b, z;
x = (((sign = (x < 0.0)) != 0) ? -x : x);
l = 0;
if (x >= 4.0) {
l = -1;
x = 1.0 / x;
y = 0;
goto atl;
} else {
if (x < 0.25) {
y = 0;
goto atl;
}
}
y = aint(x / 0.5);
z = y * 0.5;
x = (x - z) / (x * z + 1);
atl:
z = x * x;
b = ((((893025.0 * z + 49116375.0) * z + 425675250.0) * z +
1277025750.0) * z + 1550674125.0) * z + 654729075.0;
a = (((13852575.0 * z + 216602100.0) * z + 891080190.0) * z +
1332431100.0) * z + 654729075.0;
a = (a / b) * x + atanc[y];
if (l)
a=piover2 - a;
return sign ? -a : a;
}
/* atan2(y,x) Return arctangent in radians of y/x,
range -pi to pi */
double INTRIG_atan2(y, x)
double y, x;
{
double temp;
if (x == 0.0) {
if (y == 0.0) /* Special case: atan2(0,0) = 0 */
return 0.0;
else if (y > 0)
return piover2;
else
return -piover2;
}
temp = INTRIG_atan(y / x);
if (x < 0.0) {
if (y >= 0.0)
temp += pic;
else
temp -= pic;
}
return temp;
}
/* asin(x) Return arcsine in radians of x */
double INTRIG_asin(x)
double x;
{
if (INTRIG_fabs(x)>1.0) {
fprintf(stderr,
"\nInverse trig functions lose much of their gloss when");
fprintf(stderr,
"\ntheir arguments are greater than 1, such as the");
fprintf(stderr,
"\nvalue %g you passed.\n", x);
exit(1);
}
return INTRIG_atan2(x, INTRIG_sqrt(1 - x * x));
}
/* Perform ray trace in specific spectral line */
static trace_line(line, ray_h)
int line;
double ray_h;
{
int i;
object_distance = 0.0;
ray_height = ray_h;
from_index = 1.0;
for (i = 1; i <= current_surfaces; i++) {
radius_of_curvature = s[i][1];
to_index = s[i][2];
if (to_index > 1.0)
to_index = to_index + ((spectral_line[4] -
spectral_line[line]) /
(spectral_line[3] - spectral_line[6])) * ((s[i][2] - 1.0) /
s[i][3]);
transit_surface();
from_index = to_index;
if (i < current_surfaces)
object_distance = object_distance - s[i][4];
}
}
static INTRIG_trace_line(line, ray_h)
int line;
double ray_h;
{
int i;
object_distance = 0.0;
ray_height = ray_h;
from_index = 1.0;
for (i = 1; i <= current_surfaces; i++) {
radius_of_curvature = s[i][1];
to_index = s[i][2];
if (to_index > 1.0)
to_index = to_index + ((spectral_line[4] -
spectral_line[line]) /
(spectral_line[3] - spectral_line[6])) * ((s[i][2] - 1.0) /
s[i][3]);
INTRIG_transit_surface();
from_index = to_index;
if (i < current_surfaces)
object_distance = object_distance - s[i][4];
}
}
static transit_surface() {
double iang, /* Incidence angle */
rang, /* Refraction angle */
iang_sin, /* Incidence angle sin */
rang_sin, /* Refraction angle sin */
old_axis_slope_angle, sagitta;
if (paraxial) {
if (radius_of_curvature != 0.0) {
if (object_distance == 0.0) {
axis_slope_angle = 0.0;
iang_sin = ray_height / radius_of_curvature;
} else
iang_sin = ((object_distance -
radius_of_curvature) / radius_of_curvature) *
axis_slope_angle;
rang_sin = (from_index / to_index) *
iang_sin;
old_axis_slope_angle = axis_slope_angle;
axis_slope_angle = axis_slope_angle +
iang_sin - rang_sin;
if (object_distance != 0.0)
ray_height = object_distance * old_axis_slope_angle;
object_distance = ray_height / axis_slope_angle;
return;
}
object_distance = object_distance * (to_index / from_index);
axis_slope_angle = axis_slope_angle * (from_index / to_index);
return;
}
if (radius_of_curvature != 0.0) {
if (object_distance == 0.0) {
axis_slope_angle = 0.0;
iang_sin = ray_height / radius_of_curvature;
} else {
iang_sin = ((object_distance -
radius_of_curvature) / radius_of_curvature) *
sin(axis_slope_angle);
}
iang = asin(iang_sin);
rang_sin = (from_index / to_index) * iang_sin;
old_axis_slope_angle = axis_slope_angle;
axis_slope_angle = axis_slope_angle + iang - asin(rang_sin);
sagitta = sin((old_axis_slope_angle + iang) / 2.0);
sagitta = 2.0 * radius_of_curvature*sagitta*sagitta;
object_distance = ((radius_of_curvature * sin(
old_axis_slope_angle + iang)) *
cot(axis_slope_angle)) + sagitta;
return;
}
rang = -asin((from_index / to_index) *
sin(axis_slope_angle));
object_distance = object_distance * ((to_index *
cos(-rang)) / (from_index *
cos(axis_slope_angle)));
axis_slope_angle = -rang;
}
static INTRIG_transit_surface() {
double iang, /* Incidence angle */
rang, /* Refraction angle */
iang_sin, /* Incidence angle sin */
rang_sin, /* Refraction angle sin */
old_axis_slope_angle, sagitta;
if (paraxial) {
if (radius_of_curvature != 0.0) {
if (object_distance == 0.0) {
axis_slope_angle = 0.0;
iang_sin = ray_height / radius_of_curvature;
} else
iang_sin = ((object_distance -
radius_of_curvature) / radius_of_curvature) *
axis_slope_angle;
rang_sin = (from_index / to_index) * iang_sin;
old_axis_slope_angle = axis_slope_angle;
axis_slope_angle = axis_slope_angle + iang_sin - rang_sin;
if (object_distance != 0.0)
ray_height = object_distance * old_axis_slope_angle;
object_distance = ray_height / axis_slope_angle;
return;
}
object_distance = object_distance * (to_index / from_index);
axis_slope_angle = axis_slope_angle * (from_index / to_index);
return;
}
if (radius_of_curvature != 0.0) {
if (object_distance == 0.0) {
axis_slope_angle = 0.0;
iang_sin = ray_height / radius_of_curvature;
} else {
iang_sin = ((object_distance -
radius_of_curvature) / radius_of_curvature) *
INTRIG_sin(axis_slope_angle);
}
iang = INTRIG_asin(iang_sin);
rang_sin = (from_index / to_index) * iang_sin;
old_axis_slope_angle = axis_slope_angle;
axis_slope_angle = axis_slope_angle +
iang - INTRIG_asin(rang_sin);
sagitta = INTRIG_sin((old_axis_slope_angle + iang) / 2.0);
sagitta = 2.0 * radius_of_curvature*sagitta*sagitta;
object_distance = ((radius_of_curvature * INTRIG_sin(
old_axis_slope_angle + iang)) *
INTRIG_cot(axis_slope_angle)) + sagitta;
return;
}
rang = -INTRIG_asin((from_index / to_index) *
INTRIG_sin(axis_slope_angle));
object_distance = object_distance * ((to_index *
INTRIG_cos(-rang)) / (from_index *
INTRIG_cos(axis_slope_angle)));
axis_slope_angle = -rang;
}
/* Initialise when called the first time */
autodesk_init() {
int i, j;
spectral_line[1] = 7621.0; /* A */
spectral_line[2] = 6869.955; /* B */
spectral_line[3] = 6562.816; /* C */
spectral_line[4] = 5895.944; /* D */
spectral_line[5] = 5269.557; /* E */
spectral_line[6] = 4861.344; /* F */
spectral_line[7] = 4340.477; /* G'*/
spectral_line[8] = 3968.494; /* H */
/* Load test case into working array */
clear_aperture = 4.0;
current_surfaces = 4;
for (i = 0; i < current_surfaces; i++)
for (j = 0; j < 4; j++)
s[i + 1][j + 1] = testcase[i][j];
}
autodesk() {
double od_fline, od_cline;
int niter = 1000; /* Iteration counter */
/* Perform ray trace the specified number of times. */
for (itercount = 0; itercount < niter; itercount++) {
for (paraxial = 0; paraxial <= 1; paraxial++) {
/* Do main trace in D light */
trace_line(4, clear_aperture / 2.0);
od_sa[paraxial][0] = object_distance;
od_sa[paraxial][1] = axis_slope_angle;
}
paraxial = FALSE;
/* Trace marginal ray in C */
trace_line(3, clear_aperture / 2.0);
od_cline = object_distance;
/* Trace marginal ray in F */
trace_line(6, clear_aperture / 2.0);
od_fline = object_distance;
aberr_lspher = od_sa[1][0] - od_sa[0][0];
aberr_osc = 1.0 - (od_sa[1][0] * od_sa[1][1]) /
(sin(od_sa[0][1]) * od_sa[0][0]);
aberr_lchrom = od_fline - od_cline;
max_lspher = sin(od_sa[0][1]);
/* D light */
max_lspher = 0.0000926 / (max_lspher * max_lspher);
max_osc = 0.0025;
max_lchrom = max_lspher;
}
}
INTRIG_autodesk() {
double od_fline, od_cline;
int niter = 1000; /* Iteration counter */
/* Perform ray trace the specified number of times. */
for (itercount = 0; itercount < niter; itercount++) {
for (paraxial = 0; paraxial <= 1; paraxial++) {
/* Do main trace in D light */
INTRIG_trace_line(4, clear_aperture / 2.0);
od_sa[paraxial][0] = object_distance;
od_sa[paraxial][1] = axis_slope_angle;
}
paraxial = FALSE;
/* Trace marginal ray in C */
INTRIG_trace_line(3, clear_aperture / 2.0);
od_cline = object_distance;
/* Trace marginal ray in F */
INTRIG_trace_line(6, clear_aperture / 2.0);
od_fline = object_distance;
aberr_lspher = od_sa[1][0] - od_sa[0][0];
aberr_osc = 1.0 - (od_sa[1][0] * od_sa[1][1]) /
(INTRIG_sin(od_sa[0][1]) * od_sa[0][0]);
aberr_lchrom = od_fline - od_cline;
max_lspher = INTRIG_sin(od_sa[0][1]);
/* D light */
max_lspher = 0.0000926 / (max_lspher * max_lspher);
max_osc = 0.0025;
max_lchrom = max_lspher;
}
}